home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 4 / Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso / Science / RLaB / toolbox / fmin.r < prev    next >
Text File  |  1994-04-25  |  3KB  |  126 lines
  1. //---------------------------------------------------------------------------
  2. //  fmin.r:
  3. //  x = fmin(f,ax,bx,tol,print)
  4. //
  5. //  find the minimum of the function f(x) in the interval [ax,bx]
  6. //
  7. //  f = The real-values function of a single variable.
  8. //
  9. //  ax,bx = left and right endpoints of search interval
  10. //
  11. //  tol = maximum length of interval of uncertainty of minimum
  12. //        tol >= 0
  13. //
  14. //  print = an optional last agrument, if nonzero a printed tracing of
  15. //          the steps is given
  16. //
  17. //  x = the value of x approximating where f attains a minimum in [ax,bx]
  18.  
  19. //  This algorithm is from: "Computer Methods for Mathematical
  20. //  Computations", Forsythe, Malcolm, and Moler, Prentice-Hall, 1976.
  21.  
  22. //  Adapted by: Duane Hanselman, University of Maine, (207)-581-2246
  23. // initialization
  24. //---------------------------------------------------------------------------
  25.  
  26. fmin = function ( f, ax, bx, tol, _print)
  27. {
  28.   local (a, b, c, d, e, fmin_data, fu, fv, fw, fx, gs, ...
  29.          num, p, print, q, r, seps, si, step, tol1, tol2, ...
  30.          u, v, w, x, xm)
  31.  
  32.   if (!exist (_print)) { print = 0; else print = _print; }
  33.   num = 1;
  34.   seps = sqrt(eps);
  35.   c = 0.5*(3.0 - sqrt(5.0));
  36.   a = ax; b = bx;
  37.   v = a + c*(b-a);
  38.   w = v; x = v;
  39.   d = 0.0; e = 0.0;
  40.   fx = f (x);
  41.   if (print) { fmin_data = [1, x, fx] ? }
  42.   fv = fx; fw = fx;
  43.   xm = 0.5*(a+b);
  44.   tol1 = seps*abs(x) + tol/3.0;
  45.   tol2 = 2.0*tol1;
  46.  
  47.   // main loop
  48.  
  49.   while (abs (x-xm) > (tol2 - 0.5*(b-a)))
  50.   {
  51.     num = num+1;
  52.     gs = 1;
  53.  
  54.     // is parabolic fit possible
  55.  
  56.     if (abs(e) > tol1)        // yes, so fit parabola
  57.     {
  58.       gs = 0;
  59.       r = (x-w)*(fx-fv);
  60.       q = (x-v)*(fx-fw);
  61.       p = (x-v)*q-(x-w)*r;
  62.       q = 2.0*(q-r);
  63.       if (q > 0.0) { p = -p; }
  64.       q = abs (q);
  65.       r = e;
  66.       e = d;
  67.  
  68.       // is parabola acceptable
  69.  
  70.       if ((abs (p) < abs (0.5*q*r)) && (p > q*(a-x)) && (p<q*(b-x)))
  71.       {
  72.         // yes, parabolic interpolation step
  73.         d = p/q;
  74.         u = x+d;
  75.         step = "   num        x        fx       parabolic";
  76.  
  77.         // f must not be evaluated too close to ax or bx
  78.  
  79.         if (((u-a) < tol2) || ((b-u) < tol2))
  80.         {
  81.           si = sign(xm-x) + ((xm-x) == 0);
  82.           d = tol1*si;
  83.         }
  84.       else    // not acceptable, must do a golden section step
  85.         gs=1;
  86.       }
  87.     }
  88.     if (gs)    // a golden-section step is required
  89.     {
  90.       if (x >= xm) { e = a-x; else e = b-x; }
  91.       d = c*e;
  92.       step = "   num        x        fx       golden   ";
  93.     }
  94.  
  95.     // f must not be evaluated too close to x
  96.  
  97.     si = sign(d) + (d == 0);
  98.     u = x + si * max ([abs (d), tol1]);
  99.     fu = f (u);
  100.     if (print) { fmin_data = [num, u, fu] ? disp(step) ? }
  101.  
  102.     // update a, b, v, w, x, xm, tol1, tol2
  103.     if (fu <= fx)
  104.     {
  105.       if (u >= x) { a = x; else b = x; }
  106.       v = w; fv = fw;
  107.       w = x; fw = fx;
  108.       x = u; fx = fu;
  109.     else    // fu > fx
  110.       if (u < x) { a = u; else b = u; }
  111.       if ( (fu <= fw) || (w == x) )
  112.       {
  113.         v = w; fv = fw;
  114.         w = u; fw = fu;
  115.       else if ( (fu <= fv) || (v == x) || (v == w) ) {
  116.         v = u; fv = fu;
  117.       }}
  118.     }
  119.     xm = 0.5*(a+b);
  120.     tol1 = seps*abs(x) + tol/3.0;
  121.     tol2 = 2.0*tol1;
  122.   }
  123.  
  124.   return x
  125. };
  126.